import sklearn
import numpy as np
from scipy import ndimage
from matplotlib import pyplot as plt
%matplotlib inline
import matplotlib.patches as mpatches
import skimage
from skimage import io
from skimage.filter import threshold_otsu
from skimage.feature import peak_local_max
from skimage.morphology import watershed
img1 = skimage.io.imread("img/cyano1.jpg")
bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (9,9) )
sub1 = plt.imshow(img1)
arrow_args = dict(arrowstyle="->")
sub1.axes.annotate('Heterocyst'
,xy= (1600,700)
,xytext = (700,300)
,xycoords = 'data'
,textcoords = 'data' #'offset points'
,arrowprops = arrow_args
,fontsize = 14
)
sub1.axes.annotate('Vegetative'
,xy = (1080,960)
,xytext = (400,600)
,xycoords = 'data'
,textcoords = 'data' #'offset points'
,arrowprops = arrow_args
,fontsize = 14
)
def boundbox1(tryit):
# np.nonzero(np.apply_over_axes( np.argmax, thelabel,[1]))[0].min()
minX = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[1]))[0].min()
maxX = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[1]))[0].max()
minY = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[0]))[1].min()
maxY = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[0]))[1].max()
return [minX,maxX,minY,maxY]
def buildrect(arr):
minY,maxY,minX,maxX = arr
rect1 = mpatches.Rectangle(xy=(minX,minY)
,width = maxX - minX
,height = maxY - minY
,fill = False
,edgecolor = 'red'
,linewidth = 2)
return rect1
def img_to_labels1(inp_img):
thresh = threshold_otsu(inp_img)
otsu0 = inp_img < thresh
dist0 = ndimage.distance_transform_edt(otsu0)
j=34
peaks0 = peak_local_max(dist0 ,indices = False
,footprint=np.ones((j,j)))
print peaks0.cumsum()
marker0 = ndimage.label(peaks0)[0]
labels0 = watershed(-dist0, marker0, mask=otsu0)
# plt.imshow(labels0)
# plt.show()
return labels0
def make_windows(inp_labels, inp_img):
bb0 = []
for i in range(1,inp_labels.max() + 1):
label_i = (inp_labels == i)
bb0.append(boundbox1(label_i))
rect_arr = [buildrect(bb0[i]) for i in np.arange(bb0.__len__())]
fig = None
s0 = None
fig, (s0) = plt.subplots(nrows = 1, ncols = 1, figsize = (12,12))
s0.imshow(inp_img, cmap='gray')
for rect in rect_arr:
s0.add_patch(rect)
fig.show()
return fig
def algo1(inp_img):
labels = img_to_labels1(inp_img)
out_fig = make_windows(labels,inp_img)
return out_fig
img1_gray = skimage.color.rgb2gray(img1)
plt.imshow(img1_gray, cmap='gray')
cross = img1[1200:2000,0:500]
cross = skimage.color.rgb2gray(cross)
algo1(cross)
inp_img = cross
thresh = threshold_otsu(inp_img)
otsu0 = inp_img < thresh
dist0 = ndimage.distance_transform_edt(otsu0)
j=34
peaks0 = peak_local_max(dist0 ,indices = False
,footprint=np.ones((j,j)))
print peaks0.cumsum()
img = peaks0
peaks_array =[]
peak_counter = 0
for y in range(0,img.shape[0]):
for x in range(0,img.shape[1]):
if (img[y,x] == True):
peak_counter = peak_counter + 1
peaks_array.append([peak_counter,y,x])
img = cross
offset = 40
bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (25,25) )
sub1 = plt.imshow(img, cmap='gray')
arrow_args = dict(arrowstyle="->")
for i in peaks_array:
sub1.axes.annotate('Peak '+ str(i[0] - 1)
,xy= (i[2],i[1])
,xytext = (i[2]+offset,i[1]+offset)
,xycoords = 'data'
,textcoords = 'data' #'offset points'
,arrowprops = arrow_args
,fontsize = 12
)
X = np.array(peaks_array)[:,(1,2)]
from sklearn.neighbors import kneighbors_graph
weighted_graph = kneighbors_graph(X, 5, mode='distance')
def strand_v1(dist, pres, past):
"""
args:
dist: sklearn distance kneighbors_graph
pres: present label
past: past label
"""
m, n = dist.shape
n1nbr = None
for i in range(n):
if i in past or dist[pres, i]==0:
continue
if n1nbr is None or dist[pres,i]<=dist[pres,n1nbr]:
n1nbr=i
# print pres, past, n1nbr
if len(past) == 0 or n1nbr is not None:
return strand_v1(dist, n1nbr, past+[pres])
return past + [pres]
def strand_v2(dist, pres, past):
"""
Uses a mean std deviation metric to ensure neighbors are cells reasonably close
to one another, based on
args:
dist: sklearn distance kneighbors_graph
pres: present label
past: past label
"""
m, n = dist.shape
n1nbr = None
for i in range(n):
if i in past or dist[pres, i]==0:
continue
if n1nbr is None or dist[pres,i]<=dist[pres,n1nbr]:
n1nbr=i
# print pres, past, n1nbr
if len(past) == 0 or n1nbr is not None:
rawr = dist[]
if n1nbr < and np.std(past+[pres])==0:
return strand_v2(dist, n1nbr, past+[pres])
return past + [pres]
rawr = np.random.rand(10)*500
print rawr
for i in range(1, len(rawr)):
print mean_std(rawr[:i])
strand_v2(weighted_graph, 3, [])
strand_v2(weighted_graph, 19, [])
strand_v1(weighted_graph, 15, [])
strand_v1(weighted_graph, 5, [])
hand_labeled = [19, 13, 12, 10, 9, 7, 8]
# labels=tuple(strand_v1(weighted_graph, 19, []))
labels = hand_labeled
print labels
X[labels,:]
coefficients = np.polyfit(X[labels,1], X[labels,0], len(labels)-1)
print coefficients
polynomial = np.poly1d(coefficients)
ys = polynomial(X[labels,1])
fig,ax=plt.subplots(figsize=(15,15))
ax.imshow(img, cmap='gray')
ax.plot(X[labels,1], ys, 'r', lw=3)
ax.set_ybound(img.shape[0],0)
ax.set_xbound(img.shape[1],0)
img = cross
offset = 40
bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (25,25) )
sub1 = plt.imshow(img, cmap='gray')
arrow_args = dict(arrowstyle="->")
for i in peaks_array:
sub1.axes.annotate('Peak '+ str(i[0] - 1)
,xy= (i[2],i[1])
,xytext = (i[2]+offset,i[1]+offset)
,xycoords = 'data'
,textcoords = 'data' #'offset points'
,arrowprops = arrow_args
,fontsize = 12
)
# east_strand = strand_v1(weighted_graph, 19, [])
east_strand = hand_labeled
north_bound = east_strand + [11, 6]
print north_bound
west_bound = east_strand + [11, 14]
print west_bound
south_bound = east_strand + [11, 17]
print south_bound
def show_path_and_der(coords, labels, do_plot=True):
labels = tuple(labels)
coefficients = np.polyfit(coords[labels,1], coords[labels,0], len(labels)-1)
polynomial = np.poly1d(coefficients)
ys = polynomial(coords[labels,1])
dcoefficients = np.polyder(coefficients)
print dcoefficients
dpolynomial = np.poly1d(dcoefficients)
dys = dpolynomial(coords[labels,1])
# d2coefficients =np.polyder(coefficients, 2)
# print d2coefficients
# d2polynomial = np.poly1d(d2coefficients)
# d2ys = d2polynomial(coords[labels,1])
if do_plot:
fig,ax=plt.subplots(figsize=(10,10))
ax.imshow(img, cmap='gray')
ax.plot(coords[labels,1], ys, 'r', lw=3)
ax.plot(coords[labels,1], dys, 'g', lw=3)
# ax.plot(coords[labels,1], d2ys, 'b', lw=3)
# ax.set_ybound(img.shape[0],0)
# ax.set_xbound(img.shape[1],0)
show_path_and_der(X, north_bound[:-1])
show_path_and_der(X, north_bound)
show_path_and_der(X, west_bound[:-1])
show_path_and_der(X, west_bound)
show_path_and_der(X, south_bound[:-1])
show_path_and_der(X, south_bound)
def cost_of_path(coords, labels):
labels = tuple(labels)
def get_der_coef(coords, labels):
coefficients = np.polyfit(coords[labels,1], coords[labels,0], len(labels)-1)
dcoefficients = np.polyder(coefficients)
return dcoefficients
initial = get_der_coef(coords, labels[:-1])
decision = get_der_coef(coords, labels)[1:]
import scipy
return scipy.spatial.distance.euclidean(initial, decision)
cost_of_path(X, north_bound)
cost_of_path(X, west_bound)
cost_of_path(X, south_bound)